Analyzing topographic change in the USA

Check out the data

topo_change <- st_read('data/topo_change/topo_change_polygons.shp')
## Reading layer `topo_change_polygons' from data source `C:\Users\keman\Dropbox\Documents\Courses\WR674\Mining Impacts\mining_impacts_to_elevation\data\topo_change\topo_change_polygons.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8958 features and 51 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -2332814 ymin: 407783 xmax: 2238422 ymax: 3128340
## epsg (SRID):    NA
## proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
#Checkout the column names
#names(topo_change)
biggest_area <- topo_change %>%
  arrange(desc(AREA_SQ_KM)) %>%
  slice(1:10)
# Map all mines
mapview(biggest_area)

Checkout mine in Arizona

# Subset to just Arizona
az_mine <- biggest_area %>%
  filter(QUADNAME == 'esperanza_mill_AZ')
#Check that it is the right site
mapview(az_mine)

Getting our own elevation datasets

Download pre-mining DEM

#Check projection of az_mine
#st_crs(az_mine)
az_raster_before <- get_elev_raster(az_mine,z=12)
## Merging DEMs
## Reprojecting DEM to original projection
## Note: Elevation units are in meters.
## Note: The coordinate reference system is:
##  +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
#Look at the structure of the data
#str(az_raster)
#Summary of the data
#summary(az_raster)
plot(az_raster_before)

#mapview(az_raster_before)

Shrink the resolution

small_az <- aggregate(az_raster_before,fact=6)
mapview(small_az)

Take the slope of the data

az_slope <- terrain(small_az,opt='slope',unit='degrees')
mapview(az_slope)

Take the aspect of the data

az_aspect <- terrain(small_az,opt='aspect',unit='degrees')
mapview(az_aspect)

3D plotting mechanism for elevation data

az_matrix <- matrix(raster::extract(small_az,
                    raster::extent(small_az),buffer=500),
                    nrow=ncol(small_az), ncol=nrow(small_az))
az_matrix %>%
  sphere_shade(texture = 'desert') %>%
  add_water(detect_water(az_matrix),color='desert') %>%
  add_shadow(ray_shade(az_matrix,zscale=2),0.5) %>%
  add_shadow(ambient_shade(az_matrix)) %>%
  plot_3d(az_matrix,zscale=10,fov = 0,theta=135,zoom=0.75)
render_snapshot(clear = TRUE)

In Class work

For this section we are really going to focus on the mine area so first I’m going to buffer our az_mine polygon by 2km and then clip both the before and after mining rasters to this new shape

#Read in raster data
az_raster_after <- raster('data/srtm_14_06.tif')
az_mine_2km <- az_mine %>%
  st_transform(2163) %>%
  st_buffer(2000) %>%
  #Reproject into SRTM datum (WGS84)
  st_transform(st_crs(az_raster_after))
#Get extent of buffered mine area
aoi <- extent(az_mine_2km)
#Crop after to area of interest
mine_after <- crop(az_raster_after,aoi)
#Crop before to area of interest (need to also reproject)
mine_before <- crop(az_raster_before %>%
                      projectRaster(.,crs=projection(mine_after)),
                    aoi) %>%
  projectRaster(.,mine_after) #match resolution
#Force rasters to be exactly same size
mine_before_same <- trim(mine_before)
mine_after_same <- crop(mine_after,mine_before_same)
#Stack these things in a rasterBrick (which we can do because we 
#made them the exact same size)
mines <- brick(mine_before_same,mine_after_same)
#Add a dem_difference layer
mines[[3]] <- mines[[1]] - mines[[2]]
#Rename raster layers
names(mines) <- c('Pre_mining','Post_mining','DoD')
plot(mines)

All data for this section should be the data clipped to the mining area (generated in the code chunk above)

1) Get the slope of the SRTM (post-mining) dataset and plot this data (use the data that is clipped to the mine only).

az_slope_after <- terrain(mine_after_same,opt='slope',unit='degrees')
mapview(az_slope_after)

2) Use the brick code from above to plot slopes for the pre-mining, post-mining, and difference of slope

az_slope_before <- terrain(mine_before_same,opt='slope',unit='degrees')
az_slope_diff<-terrain(mines[[3]], opt='slope', unit='degrees')
slope<-brick(az_slope_after,az_slope_before,az_slope_diff)
names(slope) <- c('Post_mining','Pre_mining','DoD')
plot(slope)

3) Use rayshader to make a 3d plot of the pre-mining mine

az_before_matrix <- matrix(raster::extract(mine_before_same,
                    raster::extent(mine_before_same)),
                    nrow=ncol(mine_before_same), ncol=nrow(mine_before_same))
az_before_matrix %>%
  sphere_shade(texture = 'desert') %>%
  add_water(detect_water(az_before_matrix),color='desert') %>%
  add_shadow(ray_shade(az_before_matrix,zscale=2),0.5) %>%
  add_shadow(ambient_shade(az_before_matrix)) %>%
  plot_3d(az_before_matrix,zscale=10,fov = 0,theta=135,zoom=0.75)
render_snapshot(clear = TRUE)

4) Use rayshader to make 3d plot of the post-mining mine

az_after_matrix <- matrix(raster::extract(mine_after_same,
                    raster::extent(mine_after_same)),
                    nrow=ncol(mine_after_same), ncol=nrow(mine_after_same))

az_after_matrix %>%
  sphere_shade(texture = 'desert') %>%
  add_water(detect_water(az_after_matrix),color='desert') %>%
  add_shadow(ray_shade(az_after_matrix,zscale=2),0.5) %>%
  add_shadow(ambient_shade(az_after_matrix)) %>%
  plot_3d(az_after_matrix,zscale=10,fov = 0,theta=135,zoom=0.75)
render_snapshot(clear = TRUE)